#ifndef GEOMETRY_MEDIAN_HPP
#define GEOMETRY_MEDIAN_HPP

#include "point.hpp"
#include "line.hpp"

/* Returns the geometric median of n points, i.e. a point that minimizes the
 * total distance to the specified points. If there are more than one point 
 * that minimizes this distance, returns any one of them.
 */
#if 0
template <class T>
Point<T> median(const Point<T> points[], size_t n)
{
}
#endif

/* Returns the geometric median of one point. */
template <class T>
Point<T> median(const Point<T> (&points)[1])
{
    return points[0];
}

/* Returns the geometric median of two points. */
template <class T>
Point<T> median(const Point<T> (&points)[2])
{
    return center(points[0], points[1]); // could as well be p[0] or p[1]
}

/* Returns the geometric median of four points. */
template <class T>
Point<T> median(const Point<T> (&points)[4])
{
    /* If the points form a convex quadrilateral, the median is the crossing
     * point of its diagonals. Otherwise, the median is one of its vertices.
     * In any case, we can just compare the six candidate points.
     */
    Point<T> M;
    if (intersect(
        LineSegment<T>(points[0], points[2]),
        LineSegment<T>(points[1], points[3]), &M) == 1)
        return M;
    if (intersect(
        LineSegment<T>(points[0], points[1]),
        LineSegment<T>(points[2], points[3]), &M) == 1)
        return M;

    /* Now that the quadrilateral is either concave or degenerate, the median
     * is one of its vertices. For simplicity, we use an exhaustive search.
     */
    T min_d = -1;
    for (int i = 0; i < 4; i++)
    {
        double d = 0;
        for (int j = 0; j < 4; j++)
        {
            if (i != j)
                d += distance(points[i], points[j]);
        }
        if (min_d < 0 || d < min_d)
        {
            min_d = d;
            M = points[i];
        }
    }
    return M;
}

#endif /* GEOMETRY_MEDIAN_HPP */
